table_patterns=read.table("table_patterns.txt",header=TRUE,colClasses="character")
uniq2=c(unique(table_patterns[,2]),"total")

thres[["thr_0.01"]]=read.table("island3to4_thr_0.01")
thres[["thr_0.02"]]=read.table("island3to4_thr_0.02")
thres[["thr_0.05"]]=read.table("island3to4_thr_0.05")
thres[["thr_0.07"]]=read.table("island3to4_thr_0.07")
thres[["thr_0.1"]]=read.table("island3to4_thr_0.1")
save(thres,file="list of island3to4.RData")

#island_001=read.table("island3to4_thr_0.01")
#island_002=read.table("island3to4_thr_0.02")
#island_005=read.table("island3to4_thr_0.05")
#island_007=read.table("island3to4_thr_0.07")
#island_01=read.table("island3to4_thr_0.1")

######################################################################################################################################################
#####resampling of combinations of 3 islands from clusters with samples from 4 islands
sum_tab_subtr=NULL
aux_problm=NULL
aux=NULL
#thres=list()
#compar=NULL
matrixCl=NULL
m=NULL
t=NULL
r=NULL
h=NULL
m_out=NULL
t_out=NULL
r_out=NULL
h_out=NULL


compar=cbind(compar,uniq2[1:42])
colnames(compar)="patterns"

for(y in 2:3)
	{
	matrixCl=thres[[y]]
	all_subtr=list()
	all_nw=vector()
	aux=which(names(matrixCl)==names(thres[y]))
	uniq1=unique(matrixCl[,aux])
	patterns=matrix(data=0,nrow=length(uniq1),ncol=length(uniq2),dimnames=list(uniq1,uniq2))
	pat1=row.names(patterns)
	pat2=colnames(patterns)
	source("functions.R")

	print(Sys.time())
	for (x in uniq1)   
		{
		subAux=matrixCl[matrixCl[,aux]==x,] 
		NbIsland=length(unique(subAux[,2]))
		count_tot=0
		if(NbIsland==4)
			{
			count_3=count_3+1
			vec=c(vec,x)
			}
			m=grep("moor",subAux[,1],value=TRUE)
			t=grep("tahi",subAux[,1],value=TRUE)
			h=grep("huah",subAux[,1],value=TRUE)
			r=grep("raia",subAux[,1],value=TRUE)
			m_out=do.call("comb_3",list(t,h,r))
			h_out=do.call("comb_3",list(m,r,t))
			r_out=do.call("comb_3",list(m,h,t))
			t_out=do.call("comb_3",list(m,h,r))
			all_comb=cbind(m_out,h_out,r_out,t_out)
			for (j in 1:dim(all_comb)[2]) 
				{
				count_tot=count_tot+1
				subtr_nw=NULL # restart subtr_nw
				subtr=subset(my3_phylo4,tips.include=all_comb[,j])
				all_subtr[[paste(x,j,sep="_")]]=subtr
				subtr_phylo=as(subtr,"phylo")
				subtr_phylo$tip.label=substr(subtr_phylo$tip.label,9,9)
				subtr_phylo$edge.length=NULL
				subtr_newick=write.tree(subtr_phylo)
				subtr_nw=table_patterns[table_patterns[,1]==subtr_newick,2]
				sumary_subtr=c(paste(x,j,sep="_"),subtr_newick,subtr_nw)
				if (length(sumary_subtr)==2)
					{
					sumary_subtr=c(sumary_subtr, "pb")
					}
				sum_tab_subtr=rbind(sum_tab_subtr,sumary_subtr)
				patterns[pat1==x,pat2==subtr_nw]=patterns[pat1==x,pat2==subtr_nw]+1				}
			patterns[pat1==x,pat2=="total"]=count_tot
			}
		}
	write.table(patterns,paste("patterns for ",names(thres[y]),"_subset3of4_",Sys.Date(),".txt"), sep="\t",row.names=TRUE,col.names=TRUE)
	write.table(sum_tab_subtr,paste("summary of table subtree",names(thres[y]),"_subset3of4.txt"))
save(all_subtr,file=paste("list of all subtrees ",names(thres[y]),"_subset3of4_",Sys.Date(),".RData"))
	print(Sys.time())
	}



####funtioning and copied to patterns topology (24.11.2010)
#########################################################################################################################################################
###counting patterns for each threshold
sum_tab_subtr=NULL
aux_problm=NULL
aux=NULL
compar=NULL
thres=list()

thres[["thr_0.01"]]=read.table("island3to4_thr_0.01")
thres[["thr_0.02"]]=read.table("island3to4_thr_0.02")
thres[["thr_0.05"]]=read.table("island3to4_thr_0.05")
thres[["thr_0.07"]]=read.table("island3to4_thr_0.07")
thres[["thr_0.1"]]=read.table("island3to4_thr_0.1")
 
thres[["thr_0.01"]]=read.table("subset_test_island3to4_thr001.txt")

for(y in 1:length(thres))
{
matrixCl=thres[[y]]
all_subtr=list()
all_nw=vector()
aux=which(names(matrixCl)==names(thres[y]))
uniq1=unique(matrixCl[,aux])
patterns=matrix(data=0,nrow=length(uniq1),ncol=length(uniq2),dimnames=list(uniq1,uniq2))
pat1=row.names(patterns)


for (x in uniq1)   
{
	subAux=matrixCl[matrixCl[,aux]==x,] 
	NbIsland=length(unique(subAux[,2])) 
	aux_1=combn(as.vector(subAux[,1]),NbIsland)
	count_tot=0
 	for (j in 1:dim(aux_1)[2]) 
	{
		test=substr(aux_1[,j],9,9)
		if (length(unique(test))==NbIsland)
		{
			count_tot=count_tot+1
			subtr_nw=NULL # restart subtr_nw
			subtr=subset(my3_phylo4,tips.include=aux_1[,j])
			all_subtr[[paste(x,j,sep="_")]]=subtr
			subtr_phylo=as(subtr,"phylo")
			subtr_phylo$tip.label=substr(subtr_phylo$tip.label,9,9)
			subtr_phylo$edge.length=NULL
			subtr_newick=write.tree(subtr_phylo)
			subtr_nw=table_patterns[table_patterns[,1]==subtr_newick,2]
			sumary_subtr=c(paste(x,j,sep="_"),subtr_newick,subtr_nw)
			if (length(sumary_subtr)==2)
			{
				sumary_subtr=c(sumary_subtr, "pb")
			}
			sum_tab_subtr=rbind(sum_tab_subtr,sumary_subtr)
			patterns[pat1==x,pat2==subtr_nw]=patterns[pat1==x,pat2==subtr_nw]+1
		}			
	} 
	patterns[pat1==x,pat2=="total"]=count_tot
}
write.table(patterns,paste("patterns for ",names(thres[y]),"_subset_",Sys.Date(),".txt"), sep="\t",row.names=TRUE,col.names=TRUE)
write.table(sum_tab_subtr,paste("summary of table subtree",names(thres[y]),"_subset.txt"))
save(all_subtr,file=paste("list of all subtrees ",names(thres[y]),"_subset_",Sys.Date(),".RData"))

count_cols_thr=apply(patterns[,(1:42)],2,sum)
count_rows_thr=apply(patterns[,(1:42)],1,sum)
t=count_cols_thr
t=cbind(t,names(t))
colnames(t)=c(names(thres[y]),"patterns")
if (compar=NULL) {compar=cbinb(compar,t[,2])}
else { compar=merge(compar,t)}
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
table_comp=merge(compar,topo)
}


count_cols_thr1=apply(patterns[,(1:42)],2,sum)
count_rows_thr1=apply(patterns[,(1:42)],1,sum)
t1=count_cols_thr1
t1=cbind(t1,names(t1))
colnames(t1)=c("thr001","patterns")


count_cols_thr2=apply(patterns[,(1:42)],2,sum)
count_rows_thr2=apply(patterns[,(1:42)],1,sum)
t2=count_cols_thr2
t2=cbin(t2,names(t2))
colnames(t2)=c("thr002","patterns")
compar_t1t2=merge(t1,t2)
cmpar_t1t2=compar_t1t2[,-2]
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
table_comp=merge(compar_t1t2,topo)


count_cols_thr5=apply(patterns[,(1:42)],2,sum)
count_rows_thr5=apply(patterns[,(1:42)],1,sum)
t5=count_cols_thr5
t5=cbind(t5,names(t5))
colnames(t5)=c("thr005","patterns")
compar_t1t2t5=merge(compar_t1t2,t5)
topo=read.table("topolgy_info_or_not_01.txt",header=TRUE)
topo=topo[,-5]
table_comp=merge(compar_t1t2t5,topo)


}

####counting total events for each pattern (columns) and for each cluster (rows)
count_cols_thr1=apply(patterns[,(1:42)],2,sum)
count_rows_thr1=apply(patterns[,(1:42)],1,sum)
t1=count_cols_thr1
t1=cbind(t1,names(t1))
colnames(t1)=c("thr001","patterns")

####test if column total and sum of row are equal

for(m in names(count_rows))
	{
	if(patterns[row.names(patterns)==m,43]==count_rows[names(count_rows)==m])
	print("no problem")
		else{ print(paste("problem with",m,sep="\t"))
			aux_problm=c(aux_problm,m)}
	}
	
	
